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the  soil.  The  new  model  was  evaluated  for  its  accuracy  in  reproducing  free-field  stresses  for  clcy  and  sand, 
and  can  be  used  immediately  adjacent  to  the  explosive  source.  Due  to  large  deformations  and  very  high 
strains  and  strain  rates  in  the  soil  near  the  charge,  it  was  difficult  to  maintain  the  stability  of  the  calculation. 
A  method  of  maintaining  this  stability  without  adversely  affecting  the  character  of  the  free-field  stress  and 
velocity  time-histories  was  developed  and  validated.  A  suitable  method  for  simulating  a  nonreflecting 
boundary  for  soil  responding  inelastically  was  also  developed.  Calculations  were  performed  for  sand  and 
clay  soils  to  verify  that  DYNA3D  with  all  of  the  modifications  was  capable  of  adequately  predicting  free- 
field  stresses  and  velocities.  These  predictions  agreed  reasonably  well  with  test  data  and  indicated  that  the 
modified  version  of  DYNA3D  would  be  useful  in  investigating  SSI.  The  nonreflecting  boundary  methodol¬ 
ogy  was  also  verified. 

One-dimensional  cylindrical  and  spherical  geometry  assumptions  were  imposed  and  analyzed  to  deter¬ 
mine  if  it  was  reasonable  to  perform  the  SSI  study  using  a  two-dimensional  (2-D)  plane  strain  model.  These 
analyses  indicated  that  a  charge  can  be  selected  for  the  plane  strain  calculation  so  that  the  peak  stress,  maxi¬ 
mum  velocity,  and  stress  and  velocity  gradients  at  the  center  of  the  structure  match  those  which  are  pre¬ 
dicted  using  the  correct  charge  and  a  three-dimensional  geometry.  A  2-D  model  was  selected  to  study  SSI. 
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This  raport  documentB  the  analyses  performed  to  develop  a  methodology 
for  studying  soil-structure  interaction  when  a  conventional  weapon  detonates 
near  a  buried  hardened  structure.  The  work  was  sponsored  and  funded  through 
the  In-House  Laboratory  independent  Research  Program,  work  unit  ntimber 
A91-LD-003,  work  unit  title,  "Parameters  Affecting  Loads  on  Buried  Structures 
Subjected  to  Localized  Blast  Effects." 

These  analyses  were  performed  in  the  Structures  Laboratory  (SL),  U.S. 
Army  Engineer  Waterways  Experiment  Station  (WES),  by  Mr.  James  T.  Bay lot. 
Research  Structural  Engineer,  Structural  Mechanics  Division  (SMD),  under  the 
general  supervision  of  Messrs.  Bryant  Mather  and  James  T.  Ballard,  Director 
and  Assistant  Director,  SL,  respectively;  and  under  the  direct  supervision  of 
Drs.  Jimmy  P.  Balsara,  Chief,  SMD,  and  Robert  Hall,  Chief,  Analysis  Group, 
SMD.  This  report  was  prepared  by  Mr.  Baylot.  Mr.  Leonard  I.  Huskey  is 
manager  of  the  ILIR  Program,  and  Ms.  Mary  K.  Vincent,  Chief,  Office  of 
Technical  Programs  and  Plans,  provides  overall  coordination  for  the  program. 

At  the  time  of  publication  of  this  report.  Director  of  WES  was 
Dr.  Robert  W.  Whalin.  Commander  and  Deputy  Director  was  COL  Leonard  G. 
Hassell,  EN. 
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CONVFRSION  FACTORS,  NON-SI  TO  SI  < METRIC)  UNITS  OF  MEASUREMENT 
Non-SI  units  of  measurenant  used  in  this  report  can  be  converted  to  SI 
(metric)  units  as  follows: 


Multlplv 

BX _ 

feet 

0.3048 

inches 

25.4 

pounds 

(force) 

4.47222 

pounds 

(force)  per  square  inch 

0.006894757 

pounds 

(mass) 

0.4535924 

pounds 

(mass)  per  cubic  foot 

16.01846 

To  Obtain 

metres 

millimetres 

newtons 

megapascals 

kilograms 

kilograms  per 
ciibic  metre 
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Parameters  Affecting  Loads  on  Buried  Structurea 
Sublected  to  Localized  Blast  Effects 


PART  I :  INTRODUCTION 

Backaroiind 

1 .  The  military  has  many  important  facilities  which  must  be  capable  of 
surviving  a  conventional  weapons  attack.  Quite  often  these  facilities  are 
massive  reinforced  concrete  rectangular  structxires  in  a  buried  configuration. 
These  structures  are  extremely  expensive  to  construct,  and  good  design 
procedures  are  needed  to  provide  an  economical,  yet  safe  design.  Since  these 
facilities  will  be  located  at  sites  throughout  the  world,  there  are  many 
different  types  of  backfill  materials  which  may  be  used. 

2.  There  are  two  options  availcUale  to  those  responsible  for  designing 
these  structiires.  They  can  import  a  favorable  soil  type,  such  as  sand,  to  be 
used  as  the  backfill  material,  or  they  can  use  the  native  soil.  The  cost  of 
importing  the  sand  is  offset,  at  least  partially,  by  the  reduced  structure 
cost.  Once  both  designs  are  considered,  the  cheaper  method  can  be  selected. 
However,  since  there  is  less  confidence  in  designing  the  structure  in  the 
other  backfill  types,  the  designer  will  often  select  importing  the  sand 
backfill  material.  This  leads  to  structures  that  are  much  more  expensive  to 
build  than  they  would  have  been  if  more  reliable  methods  of  assessing  the 
importance  of  soil  porameters  were  available. 

3.  Probably  the  most  difficult  phase  of  designing  a  buried  structure  to 
resist  the  effects  of  a  conventional  weapons  detonation  is  the  deteinnination 
of  the  loads  on  the  structure.  Generally  the  free-field  stresses  (the 
stresses  which  would  be  present  if  the  structure  was  absent)  at  the  location 
of  the  point  on  the  structure  are  computed.  These  stresses  are  modified  to 
approximate  the  effects  of  the  structure  and  its  response.  These  modified 
stresses  (interface  stresses)  are  applied  as  the  structure  loading.  Army 
Technical  Manual  (TH)  5-855-1  (U.S.  Army  Engineer  Waterways  Experiment  Station 
1986)  provides  one  method  of  determining  the  loading  from  the  free-f:leld 
stresses .  Another  method  has  been  proposed  by  Drake  et  al .  (1 987 ) . 
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4.  In  TM  5-855-1,  semiempirical  methods  to  determine  interface  stresses 
from  the  free- field  stresses  are  used.  The  method  recommended  by  Drake  uses 
continuity  of  velocities  and  stresses  at  the  interface,  along  with 
linear-elastic,  plane-wave  theory  to  predict  interface  stresses.  In  this 
method  the  loading  on  the  structure  is  predicted  using  the  free-field  stresses 
and  velocities  and  the  structxire  velocity.  This  method  will  be  referred  to  as 
the  structure-medium  interaction  (SMI)  method.  One  more  simplifying 
assumption  which  can  be  made  is  that  the  free-field  stress  is  the  acoustic 
impedance  multiplied  by  the  free-field  velocity.  Using  this  assumption,  the 
free-field  velocity  time-history  is  not  needed.  This  method  will  be  referred 
to  as  the  modified  SMI  (MSHI)  method. 

5.  A  series  of  tests  (Bay lot  et  al.  1S85)  has  been  conducted  to  study 
the  parameters  affecting  the  response  of  buried  structures  to  conventional 
weapons  effects.  These  tests  were  sponsored  jointly  by  the  Air  Force 
Engineering  and  Services  Center  (AFESC)  and  the  Office,  Chief  of  Engineers, 
U.S.  Army  (OCE)  and  will  be  referred  to  as  the  AFESC  tests.  In  the  AFESC 
tests,  eleven  tests  were  conducted  in  a  sauid  backfill  material.  The  tests 
were  designed  to  study  the  effects  of  span- to- thickness  ratio  and 
reinforcing-steel  ratio  of  the  slab  and  charge  orientation  and  staxidoff 
distance . 

€.  Data  from  these  tests  were  studied  (Bay lot  and  Hayes  1989)  in  great 
detail  to  evaluate  the  different  methods  of  predicting  interface  stresses  from 
free-field  stresses  and  velocities.  The  methods  presented  in  TM  5-855-1  of 
predicting  interface-pressure  loads  from  free-field  stresses  significantly 
overpredict  the  loading  on  the  structure  and,  thus,  structural  response  is 
significantly  overpredicted.  The  SMI  and  MSMI  methods  appear  to  do  a 
reasonable  job  of  predicting  the  very  early  time-interface  stresses,  but  may 
significantly  underpredict  later  time  loading.  This  later  time  loading  can 
contribute  significantly  to  the  response  of  the  structure,  especially  when 
deformations  are  large.  For  thinner  slabs,  these  data  indicate  that  some 
other  method  of  predicting  the  later  time  loading  is  needed. 

7.  The  Defense  Nuclear  Agency  recently  sponsored  a  series  of 
conventional  weapons  backfill  experiments  (CONWEB)  to  study  the  effect  of  soil 
type  on  structure  response  (Hayes  1989).  They  showed  that  the  soil  type  is 
extremely  important  in  determining  the  loads  on  and,  thus,  the  response  of 


5 


buried  structures.  In  a  clay  backfill  experiment,  the  structure  failed;  while 
in  an  identical  experiment  in  a  sand  backfill,  minor  damage  was  incurred. 

8.  These  recent  studies  indicate  that  one  current  method  of  designing 
buried  structures  for  conventional  weapons  threats  may  be  overly  conservative, 
while  the  other  method  may  underpredict  response.  Neither  method  adequately 
predicts  the  differences  in  response  caused  by  different  backfill  materials. 
Since  the  objective  of  the  design  is  to  provide  the  most  economical  design, 
which  will  provide  the  desired  level  of  protection,  neither  method  is 
completely  satisfactory.  For  this  reason,  a  study  is  needed  to  determine 
those  parameters  which  significantly  affect  the  response  of  buried  structures. 

9.  Although  much  can  be  learned  by  conducting  experiments  on  model 
structures  and  carefully  examining  the  data,  there  are  generally  factors  which 
complicate  the  interpretation  of  results  and  make  the  type  of  study  needed 
impossible  to  perform  with  experimental  data  alone.  One  good  approach  for 
conducting  this  study  is  to  perform  calculations  to  model  the  experiments,  and 
once  these  calculations  are  performed  acceptably,  this  approach  can  be  used  to 
conduct  parameter  studies.  The  finite  element  (PE)  method  is  an  excellent 
tool  for  performing  these  analyses.  The  charge,  soil,  and  structure  can  be 
modeled  using  the  FE  method  so  that  assumptions  about  the  SSI  are  not  needed. 
After  the  calculations  have  been  performed,  the  output  may  be  studied  to  gain 
information  on  the  characteristics  of  this  SSI  and  the  parameters  that  affect 
it.  Both  two-  (2-D)  and  three-  (3-D)  dimensional  calculations  can  be  used. 

10.  FE  calculations  (Weidlinger  and  Hinman  1987)  have  been  used  to 
analyze  buried  structures  subjected  to  conventional  weapons  effects.  These 
analyses  used  equations  from  TM5-855-1  and  procedures  similar  to  those 
recommended  by  Drake,  to  determine  the  loads  on  the  structure.  Bogosian* 
used  a  "soil  island"  approach  to  analyze  the  wall  of  a  buried  structure.  In 
this  analysis,  a  small  portion  of  the  soil  in  front  of  the  structure  was 
modeled  and  stresses  were  input  on  the  free  boundary  of  the  soil.  The  "soil 
island"  method  will  be  discussed  more  fully  in  a  later  section.  Either  of 


m 

Presentation,  19  September  1989,  Mr.  David  Bogosian,  Karagozian  &  Case 
Structural  Engineers,  South  Pasadena,  CA,  at  "Conventional  Weapons  Backfill 
Test"  meeting  sponsored  by  the  Defense  Nuclear  Agency,  Alexandria,  VA. 
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these  methods  uses  assumptions  which  affect  the  loads  transferred  to  the 
structure  and  are  therefore  not  appropriate  for  the  SSI  study . 

11.  In  order  to  perform  the  SSI  analyses,  the  FE  code  must  contain 
methods  of  modeling  the  explosives  and  the  nonlinear  behavior  of  the  soil  and 
the  structxare.  The  code  should  be  capable  of  modeling  the  interface  between 
the  soil  and  structure  and  should  have  a  nonreflecting  boundary  capability. 

The  computer  code,  DYNA3D  (Hallquist  and  Benson  1987),  has  these  capabilities. 
DYNA3D  is  widely  available  at  no  cost  to  those  interested  in  performing  these 
type  of  analyses.  This,  coupled  with  the  fact  that  this  type  of  analysis  has 
not  been  performed,  indicates  that  using  DyNA3D  will  probably  be  difficult,  if 
not  impossible.  Before  DyNA3D  can  be  used  to  study  SSI,  it  must  be  shown  that 
DYNA3D  can  be  used  to  successfully  predict  free-field,  stress  and  velocity 
time-histories . 


Oblective 

12.  The  objective  of  this  study  was  to  develop  a  computational 
procedure  in  which  the  detonation  of  the  explosive,  propagation  of  stresses 
through  the  soil,  interaction  of  the  soil  with  the  structure,  and  structure 
z espouse  are  modeled  in  a  single  analysis.  Thus,  assumptions  are  not  needed 
to  compute  structure  loads  from  free-field  stresses.  This  procedure  can  then 
be  used  to  study  SSI. 


Approach 

13.  Based  on  documentation  in  its  user's  manual,  DYNA3D  should  be 
capable  of  performing  the  desired  analyses  if  the  Cap  model  used  to  model  the 
soil  can  be  modified  to  perform  in  the  very  high  stress  region  adjacent  to  the 
explosive  source.  Before  these  modifications  were  made,  calculations  were 
performed  to  determine  if  the  Cap  model  functioned  well  in  regions  of  low 
stress  for  the  clay  material  properties.  This  Cap  model  would  not  converge  to 
a  solution  for  the  sample  test  case  used,  therefore,  this  Cap  model  was  not 
suitable  for  modifications  to  make  it  perform  well  in  the  high  stress  regions. 
Pelessone  ( 1 989 )  had  also  determined  that  the  Cap  model  in  DyNA3D  was  not 
acceptable  for  calculations  in  soil  materials  and  developed  another  version  of 
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the  Cap  model  which  was  acceptable.  This  Cap  model  was  obtained,  installed  in 
DYNASD,  and  modified  so  that  it  could  be  used  in  the  very  high  stress  region 
near  the  charge. 

14.  Once  an  acceptable  Cap  model  was  developed,  DYNA3D  was  used  to 
perform  a  series  of  calculations  to  determine  if  the  free-fie.'.d  stresses  and 
velocities  could  be  adequately  predicted  for  the  CONWEB  clay  and  sand 
backfill.  One-dimensional  (1-D)  spherical  calculations  were  performed  by 
enforcing  spherical  boundary  conditions.  These  calculations  were  unstable  and 
a  method  of  stabilizing  the  calculations  without  adversely  affecting  the 
predicted  results  was  developed. 

15.  There  are  several  parameters  which  must  be  determined  in  developing 
the  FE  model  to  be  used  in  the  SSI  study.  These  parameters  could  affect  the 
accuracy  and  stability  of  the  solution  as  well  as  the  time  required  to  perform 
the  calculations.  Therefore,  1-D  spherical  calculations  were  performed  to 
determine  optimum  values  for  parameters  such  as  grid  spacing,  maximum  strain 
increment,  and  critical  time-step  ratio. 

16.  The  SSI  investigation  may  be  performed  using  a  2-D  model  of  the 
explosive,  soil,  and  structure.  As  far  as  the  charge  and  free-field  are 
concerned,  this  is  a  1-0  cylindrical  geometry  as  opposed  to  the  1-D  spherical 
geometry  of  the  free-field  if  a  3-D  model  of  the  problem  is  used.  Therefore, 
it  was  necessary  to  determine  the  effects  of  using  cylindrical  versus 
spherical  boundary  conditions. 

17.  Original  calculations  were  performed  with  the  boundary  far  enough 
away  so  that  reflections  did  not  come  back  to  the  points  of  interest  during 
the  time  of  interest.  During  the  actual  analyses  performed  to  study  SSI,  a 
large  number  of  nodes  and  elements  will  be  needed  to  adequately  model  the 
explosive,  soil,  and  structure.  For  these  computations  to  be  performed  in  a 
reasonable  amount  of  time,  the  total  number  of  nodes  and  elements  must  be 
limited.  Since  the  coarseness  of  the  grid  affects  the  calculated  results,  it 
is  desirable  to  reduce  the  number  of  nodes  and  elements  by  reducing  the  amount 
of  soil  being  modeled.  Therefore,  it  will  be  impossible  to  have  boundaries 
that  are  far  away,  and  a  nonreflecting  boundary  must  be  used  at  the  boundary. 
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Calculations  were  performed  to  assess  the  nonreflecting  boundary  used  by 
DYNA3D.  An  error  was  found  in  this  nonreflecting  boundary.  This  error  was 
corrected^  however,  the  bovindary  still  did  not  perform  adequately. 
Modifications  to  the  nonreflecting  boundary  were  made. 
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PART  XI:  MODELING  CONSIDERATIONS 

gflnargJ. 

18.  There  are  a  number  of  important  factors  which  must  be  considered  if 
a  useful  SSI  study  is  to  be  conducted.  Typically/  the  stimcture  will  be 
designed  for  a  relatively  close-in  detonation  of  the  weapon.  Since  the  weapon 
will  be  close  in,  damage  will  be  localized  and  a  higher  level  of  damage  is 
usually  acceptable.  Thus,  the  analysis  must  consider  the  closeness  of  the 
charge  to  the  structure  and  the  possibility  of  both  geometric  and  material 
nonlinearities  in  the  soil  and  structure  responses.  DYNA3D  allows  for  both 
material  and  geometric  nonlinearities. 

Explosive  Charge 

19.  The  detonation  of  the  charge  must  be  included  in  the  calculation. 
This  can  be  done  by  directly  including  the  detonation  of  the  charge  or  by 
using  a  "soil  island",  where  a  portion  of  the  soil  is  included  in  the 
calculation,  and  stresses  and/or  velocities  are  input  at  the  free  boundary  of 
the  soil.  Thus,  experimental  data  or  results  of  previous  calculations  can  be 
used  to  prescribe  the  stresses  and/or  velocities  at  the  boundary.  This  free 
soil  boundary  must  be  taken  between  the  charge  and  the  structure. 

20.  This  method  is  adequate  as  long  as  the  time  of  interest  in  the 
calculation  is  small.  Since  the  charge  is  close  to  the  structure,  the 
boundary  along  which  stresses  are  input  must  also  be  close  to  the  structure. 
When  the  stress  wave  propagating  through  the  soil  strikes  the  structure,  a 
stress  wave  is  reflected  back  into  the  soil.  This  wave  will  propagate  back 
towards  the  boundary  and  will  try  to  interact  with  the  stresses  2uid  velocities 
at  the  boundary.  Since  these  stresses  and/or  velocities  are  already 
specified,  the  interaction  cannot  occur,  and  artificial  reflections  will  occur 
off  of  the  boundary.  The  analysis  is  only  valid  until  these  reflections  reach 
the  structure  again.  Thus,  it  is  desirable  to  directly  include  the  charge  in 
the  calculation. 

21 .  DYNA3D  contains  the  Jones  Wilkins  Lee  ( JWL)  equation  of  state  model 
which  can  be  used  to  model  the  explosive  (Dobratz  1981).  DYNA3D  is  capable  of 
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modeling  the  propagation  of  the  detonation  tlirough  the  explosive  source,  with 
the  detonation  starting  from  one  or  several  points  (lines  or  planes),  or  the 
calculation  can  be  started  after  the  explosive  is  completely  detonated,  h 
study  has  been  conducted  to  determine  the  differences  between  the  two.  This 
study  showed  that  there  was  little  difference,  except  for  very  close  to 
charge,  which  is  not  of  interest  in  this  study.  Therefore,  this  study  was 
performed  neglecting  the  effects  of  the  propagation  of  the  explosion  through 
the  source.  Parameters  input  for  this  model  are  available  (Dobratx  1981)  for 
a  numl>er  of  different  explosives  types.  The  OWL  parameters  for  C-4  as  used  in 
this  study  are  listed  below.  The  mass  density  is  also  listed  laelow.  The  JWL 
equation  of  state  is:  y 

-RiV  I  -luv 

P  =  A[1-w/(R^V)le  ♦  B[1-w/(R^')]e  +  wE/V 


where  : 

P  =  Pressure,  psi 
A  =  Pressure  coefficient 
B  a  Pressure  coefficient 
=  Coefficient 
Rj  =  Coefficient 
w  a  Coefficient 
V  a  Specific  volume,  ratio  of 
original  volume 
E  a  Internal  energy 
Eg  a  Initial  internal  energy 
a  Original  mass  density 


8.044  X  10^  psi** 

1 .870  X  10®  psi 
4.5 
1 .4 
0.25 

current -to- 

.  .  3/.  3 

psi-in.  /in. 

1.305  X  10®  psi  in.Vin.^ 
1.497  X  10"^  lb  secVin/ 


* 

Letter  dated  11  May  1989,  subject:  "Analysis  of  Bum/No  Bum  Options 
for  1-D  High-Explosive  Spherical  Source  Calculation."  From  H.  D.  Zimmerman, 
California  Research  and  Technology,  Inc.,  Chattsworth,  California,  to  Dr.  j. 
G.  Jac)ison,  Chief,  (^eomechanics  Division,  Structures  Laboratory,  U.S.  Army 
Engineer  Waterways  Experiment  Station,  Vic)u5buxg,  Mississippi. 

A  table  of  factors  for  converting  non-si  units  of  measurement  to  SI 
(metric)  units  is  presented  on  page  3. 
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22.  The  nonlinear  behavior  of  the  soil  between  the  charge  auid  the 
structure  nust  be  modeled  correctly.  The  Cap  model  is  available  in  DYNA3D 
and  is  very  suitable  for  modeling  the  nonlinear  behavior  of  soils  (Sandler  and 
Rubin  1979)  and  (Simo  et  al.  1985).  Therefore,  it  was  selected  to  model  both 
the  clay  and  the  seund  experiments  in  the  CONWEB  test  series.  It  was 
determined  that  this  Cap  model  was  not  functioning  correctly,  and  another  Cap 
model  was  obtained  (Pelessone  1989).  The  new  Cap  model  is  very  similar  to  the 
one  Installed  in  DYNA/JD.  It  was  modified  so  that  the  model  could  be  used  in 
the  very  high  pressure  region  adjacent  to  the  explosive.  The  material 
properties  for  the  Cap  model  were  obtained  from  static  uniaxial  strain  and 
triaxial  compression  test  data;  however,  minor  modifications  were  needed  to 
make  them  perform  correctly.  These  modifications  will  be  discussed  later  in 
this  report.  Numerous  parametric  calculations  were  then  made  to  determine  the 
critical  element  size,  critical  time-step  size,  and  artificial  viscosity 
coefficients. 

23.  A  complete  description  of  this  Cap  model  is  provided  in  Pelessone 
(1989).  A  very  brief  summary  will  be  presented  here.  This  model  uses  soil 
mechanics  sign  conventions  (compressive  stresses  and  strains  are  positive). 

The  Cap  model  is  a  two- invariant  model  where  yielding  is  based  on  the  square 
root  of  the  second  invariant,  J2' ,  of  the  deviatoric  stress  tensor,  given  by: 

Jj'  .  1/2  s„^I 

where  s^j  is  the  deviatoric  stress  tensor  s^j  «  0^^  -  P 
o^j  is  the  stress  tensor 
P  is  the  pressure,  P  s  1^3.0 

I.,  is  the  first  invariant  of  the  stress  tensor  I.,  -  0^^ 
and  repeated  subscripts  imply  s\immation. 

J2'  is  also  given  by: 

jg'  »  1/6((0  „-0  22)^(0  22-0  33-® 

24.  The  yield  surface  is  defined  by  the  curve  in  Figure  1.  Yielding 
occurs  when y reaches  this  envelope.  This  figure  shows  that  the  yield 
surface  consists  of  thrcte  parts  designated  the  tensile  cutoff,  the  failure 
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surface  and  the  Cap.  The  tensile  cutoff  is  a  constant,  and  the  failure 
surface  is  a  function  of  and  is  given  by: 

f  s  e  -  Y  exp(-  B  i^)  +01^ 
where  8,  y,  B,  and  6  are  material  constants 
f  is  the  failure  surface 

25.  The  Cap  is  elliptical  in  shape  and  is  movable.  The  intersections 
of  the  Cap  with  the  failure  surface  and  the  axis,  respectively,  occur  at  L 
and  X,  as  shown  in  Figure  1 .  X  is  determintxl  from  the  volumetric  plastic 
strain  and  the  hardening  function  given  below.  L  is  determined  from  X  and  the 
aspect  ratio,  R,  of  the  Cap.  R  is  a  material  constant  to  be  determined  from 
test  data.  One  difference  between  the  new  Cap  and  the  one  installed  in  DYNA3D 
is  that  a  continuous  function  is  used  to  describe  the  yield  surface.  Thus, 
there  is  no  discontinuity  in  slope  at  L.  The  hardening  function  is  shown  in 
Figure  2  and  is  defined  by: 

e/  “  W[1-exp(-D1  (X-X0)-D2(X-X0)^)], 
where  e^^  is  the  volumetric  plastic  strain. 

W  is  the  maximum  volumetric  plastic  strain. 

X  is  the  current  Cap  location. 

XO  is  the  Initial  location  of  the  Cap. 

W,  D1 ,  D2,  and  XO  are  material  constants,  which  must  be  determined.  This  is 
slightly  different  from  the  Cap  installed  in  DYNA3D,  which  did  not  include  the 
term  including  the  constant,  D2. 

26.  Associative  flow  rules  are  used  to  define  the  incremental  plastic 
strains.  Thus,  the  incremental  plastic  strains  are  normal  to  the  failure 
surface.  In  this  version  of  the  Cap  model,  the  elastic  response  of  the 
material  is  defined  by  constant  bulk,  K,  and  shear,  G,  moduli.  At  very  high 
pressures,  the  bulk  modulus  increases  significantly.  For  the  material  very 
close  to  the  charge  this  increase  must  be  modeled.  The  Cap  model  was  modified 
to  include  the  following  function  for  the  elastic  bulk  modulus: 

K  *  Kj  ♦  K,  p’'2 
where  K  is  the  bulk  modulus 

Kj,  K.|,  and  are  material  constants. 
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27.  In  DYNA3D  the  function  of  the  constitutive  model  is  to  return  the 
next  stress  state  given  the  current  stress  state  and  a  strain  increment.  The 
stress-strain  response  during  this  strain  increment  may  be  highly  nonlinear. 
Some  procediure  must  be  used  to  ensure  that  stresses  outside  the  yield  surface 
are  not  predicted,  and  that  plastic  flow  is  normal  to  the  yield  surface.  In 
the  original  Cap  model,  an  iterative  procedure  was  used.  In  the  new  Cap  model 
the  strain  increment  was  subdivided  into  a  large  number  of  smaller  increments 
within  the  constitutive  model.  The  critical  strain  increment  is  given  by: 

deltep  B  0.05  (a  -  y)  [minimum  of  (1./G  and  R/9K)] 

For  the  clay  material  tested  in  the  CONWEB  series,  this  equation  implies  that 
a  strain  Increment  of  1.39  X  10~^  must  be  used.  Since  very  large  strains  are 
expected  in  the  soil  near  the  charge,  an  extremely  large  number  of  strain 
increments  will  be  needed.  For  a  large  grid  such  as  the  one  to  be  used  for 
the  SSI  analyses,  using  a  maximum  strain  increment  of  this  size  would 
drastically  increase  run  times.  Therefore,  the  Cap  model  was  changed  so  that 
the  maximum  strain  increment  is  now  a  user  input.  The  strain  increment  that 
will  be  needed  will,  in  general,  depend  on  the  material  properties  selected 
and  the  stress  state  to  which  the  material  is  subjected.  Parametric  analyses 
were  performed  to  determine  the  optimum  maximum  strain  increment  for  the  sand 
and  clay  materials. 


Artificial  Viscosity  Coefficients 

28.  Two  other  inputs  into  the  FE  model  are  the  viscosity  coefficients. 
Linear  and  quadratic  viscosity  coefficients  are  input  into  the  code  to 
stabilize  the  calculation.  In  most  FE  calculations,  the  default  values  of 
these  coefficients  are  appropriate.  They  are  large  enough  to  stabililze  the 
calculation  but  are  small  enough  so  they  do  not  adversely  affect  the  results 
of  the  calculation.  These  coefficients  are  transparent  to  the  typical  FE  code 
user  who  probably  does  not  know  that  they  are  there.  The  viscosities  are  used 
to  generate  artificial  forces  that  are  proportional  to  the  volumetric  strain 
rate  and  the  square  of  the  volumetric  strain  rate,  respectively.  The 
quadratic  viscosity  is  only  active  during  the  original  loading  of  the 
material.  These  coefficients  affect  the  character  of  the  results,  as  well  as 
the  stability  of  the  solution.  In  general,  the  smallest  amount  of  viscosity 
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which  will  stabilize  the  solution  should  be  used.  Analyses  are  needed  to 
determine  the  appropriate  viscosity  coefficients  to  be  used  in  the  SSI  study. 


Grid  Size 

29.  The  solution  will  also  be  affected  by  grid  parameters  such  as 
element  dimensions  and  location  and  types  of  boundaries  used.  Test  data  have 
shown  that  the  rise  times  associated  with  stresses  at  a  point  in  the  soil  near 
a  conventional  weapons  detonation  are  very  small,  much  smaller  than  the 
natural  period  of  the  structure  wall.  Since  this  rise  time  could  affect  the 
response  of  the  structure,  it  is  important  that  the  rise  time  is  as  close  to 
correct  as  is  possible.  In  the  calculation,  the  rise  time  is  affected  by  the 
material  properties  as  well  as  the  grid  spacing.  The  SSI  analyses  require 
that  a  large  area  be  modeled.  Thus,  the  grid  si>acing  should  be  optimized. 
Analyses  were  performed  to  determine  the  maximum  grid  spacing,  which 
adequately  predicted  the  rise  times  on  free-field  stresses. 

30.  In  the  test  events,  the  soil  around  the  structure  continued  for  a 
great  distance.  This  amount  of  soil  cannot  be  modeled  in  the  PE  analysis.  The 
FB  grid  must  be  stopped,  and  some  type  of  boundary  conditions  enforced  at  the 
boundary.  There  are  several  types  of  boundaries  that  have  been  used  to  model 
a  nonreflecting  bovindary  (Lysmer  and  Kuhlemeyer  1969)  and  (Underwood  and  Geers 
1978).  This  boundary  simulates  an  infinite  amo\int  of  soil  placed  behind  it. 
This  boundary  has  been  shown  to  be  very  effective  in  modeling  material 
responding  in  the  elastic  range  (Underwood  and  Geers  1976)  but  did  not  perform 
well  when  modeling  the  nonlinear  behavior  of  soils  (Underwood  and  Geers  1979). 
A  nonreflecting  boundary  is  available  in  DYNA3D,  and  this  boundary  could  be 
very  useful  if  it  could  be  shown  to  be  effective.  Therefore,  analyses  were 
performed  to  assess  this  nonreflecting  boundary. 

Time  Step 

31 .  Since  DYNA3D  is  an  explicit  FE  code,  the  time  step  selected  is 
critical  to  the  calculation.  DYNA3D  selects  a  critical  time  step  for  each 
element  based  on  the  dimensions  and  wave  speed  in  that  element.  The  critical 
time  step  is  the  smallest  time  step  determined  for  any  of  the  elements.  The 
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wave  speed  is  computed  based  on  constant  elastic  material  properties  and  does 
not  consider  that  these  properties  may  change  with  stress  level.  Therefore, 
it  is  possible  that  the  critical  time  step  selected  is  not  small  enough.  A 
scale  factor  to  reduce  the  critical  time  step  can  be  input  into  DYNA3D. 
Analyses  were  performed  to  determine  the  optimum  critical  time-step  factor. 

Two-Dimensional  Effects 

32.  In  the  test  event,  a  cylindrical  charge  was  placed  relatively  close 
to  the  test  structure.  The  27-in. -long  charge  was  placed  60  in.  from  the 
structure.  It  has  been  shown  that  the  free-field  stresses  and  velocities  in 
clay  and  sand,  respectively,  can  be  predicted  reasonably  well  using  a 
spherical  charge  and  assuming  spherical  symmetry  (Zimmerman  et  al.  June  1990) 
and  (Zimmerman  et  al.  October  1990).  Many  of  the  SSI  analyses  will  be 
performed  on  a  2-D  model.  In  effect,  the  charge  will  be  an  infinitely  long 
cylinder.  The  stresses  in  this  case  will  attenuate  in  a  2-D  space  rather  than 
a  3-D  space,  and  the  effects  of  this  difference  should  be  investigated. 
Axialyses  using  1-D  cylindrical  and  spherical  boundary  conditions  were 
performed  in  the  sand  and  clay  materials  to  assess  the  importance  of  this 
difference. 
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PART  III:  PARAMETER  STUDIES 


Material  Models 

33.  Constants  for  the  modified  Cap  model  have  been  determined  for  the 
CONWEB  clay*  and  sand.**  These  constants  are  based  on  static  tests  of  those 
materials.  These  constants  (Table  '• )  were  obtained  by  using  a  constitutive 
model  driver  to  simulate  static  tests  which  had  been  conducted  on  the  two 
materials . 

34.  Since  the  Cap  parameters  are  determined  using  static  tests  of 
recompacted  material,  some  modifications  may  be  needed  to  make  the  response 
match  that  of  the  material  tested  in  the  model  tests.  Nelson^  determined 
that  the  material  in  the  CONWEB  clay  test  was  stiffer  than  was  indicated  by 
the  laboratory  tests,  and  found  that  if  the  hardening  parameter,  D1 ,  was 
decreased  to  0.0004/psl,  the  calculations  would  do  a  better  job  of  matching 
data.  Therefore  in  the  clay  calculations,  a  D1  value  of  0.0004  was  used 
instead  of  that  determined  from  laboratory  tests.  The  material  properties 
using  this  value  of  D1  will  be  referred  to  as  the  modified  material  properties 
for  clay.  A  comparison  of  the  volumetric  stress-strain  cux*ves  using  the  two 
different  values  of  D1  is  given  in  Figure  3.  All  other  parameters  were  as 
determined  in  the  laboratory  tests.  No  previous  analyses  were  available  for 
the  sand  material;  therefore,  the  appropriate  material  constants  had  to  be 
determined  by  performing  analyses  aiid  comparing  them  to  test  data. 


Personal  Communication,  1  November  19B9,  Dr.  Jon  Windham,  Research  Civil 
Engineer,  U.S.  Army  Engineer  Waterways  Experiment  Station,  Vicksburg, 
Mississippi . 

mm 

Personal  Communication,  14  August  1990,  Mr.  Steve  Akers,  Research  Civil 
Engineer,  U.S.  Army  Engineer  Waterways  Experiment  Station,  Vicksburg, 
Mississippi . 

^  Letter,  13  March  1990,  from  Dr.  Ivan  Nelson,  Weidlinger  Associates,  Inc., 
Consulting  Engineers,  New  York,  New  York,  to  Dr.  Jon  Windham,  Research  Civil 
Engineer,  U.S.  Army  Engineer  Waterways  Experiment  Station,  Vicksburg, 
Mississippi . 
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Spherical  Clav  Analyses 


35.  Analyses  to  determine  optimum  values  of  the  FE  parameters  for  the 
SSI  calculation  can  be  performed  much  more  efficier.tly  using  a  1-D  spherical 
grid  as  opposed  to  using  a  full  3-0  grid,  which  would  give  the  same  results. 
Therefore  the  grid  shown  in  Figure  4  was  used.  This  grid  was  one  element 
thick  in  the  y  and  z  directions,  and  the  x  direction  represented  the  radial 
direction  of  the  spherical  calculation.  In  order  to  simulate  the  spherical 
r'jns.  the  y  and  z  dimensions  of  the  elements  must  increase  in  proportion  to 
the  X  distance  of  the  element  from  the  origin.  In  most  of  the  calculations  a 
6-degree  sector  of  a  sphere  was  used.  Boundary  conditions  were  specified  such 
that  the  nodes  must  slide  along  radial  lines  passing  through  the  origin,  since 
this  must  be  the  case  for  spherical  symmetry.  Each  of  the  elements 
representing  the  soil  was  an  8-node  solid  element.  Since  the  charge  was 
located  at  the  origin,  the  four  nodes  located  at  a  radius  of  0.0  in.  collapsed 
to  one  node,  leaving  a  five  node  solid.  The  DyNA3D  user's  manual  (Hallquist 
and  Benson  1987)  shows  that  only  4-,  6-,  and  8-node  elements  are  available  to 
model  solids;  therefore,  the  element  was  divided  into  two  4-node  elements,  as 
shown  in  Figure  5. 

36.  It  was  desired  to  compare  the  results  of  these  analyses  to  measured 
test  data  to  determine  if  the  FE  code  predicts  free-field  stresses  and 
velocities  that  are  realistic.  Test  measurements  were  made  at  3,  4,  5,  6,  and 
7  ft  away  from  the  charge,  and  20  msec  was  selected  as  the  time  of  interest. 
Therefore,  the  grid  was  made  large  enough  so  that  reflections  from  artificial 
boundaries,  such  as  the  end  of  the  grid,  did  not  reach  the  test-gage  locations 
within  20  msec.  In  most  of  the  analyses  performed,  the  soil  boundary  was 
placed  at  40  ft  from  the  center  of  the  charge. 

37.  In  the  CONWEB  experiments,  a  controlled  backfill  was  placed  in  a 
finite  test  pit  around  the  structure.  The  charge  was  buried  to  a  depth  of  5 
ft  below  the  ground  surface.  The  test  pit  boundaries  and  the  free  surface  at 
the  top  of  the  backfill  affected  the  stresses  and  velocities  measured  in  these 
experiments.  These  effects  were  not  considered  in  the  analyses,  and  this 
affects  comparisons  with  the  data. 
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38.  Typically  In  performing  FE  calculations,  an  attempt  is  made  to 
develop  a  grid  containing  elements  with  aspect  ratios  close  to  1.0.  Since  the 
y  and  z  dimensions  of  this  grid  are  growing  with  the  distance  of  the  element 
from  the  origin,  the  x  dimension  of  the  element  must  also  grow  with  distance 
from  the  origin.  In  order  for  the  first  soil  element  to  have  aspect  ratios  of 
1.0,  the  X  dimension  of  the  first  element  must  be  approximately  the  charge 
radius,  4  in.  times  sin(6),  0.4  in.  In  order  to  maintain  this  aspect  ratio, 
the  X  dimension  of  the  elements  must  grow  by  a  factor  of  1‘«-sin(6).  Using  this 
size  element  near  the  charge  and  this  growth  rate  will  produce  an  element  at 
the  5-ft  range,  which  has  an  x  dimension  of  approximately  6  in. 

39.  If  a  2-degree  sector  was  used  instead  of  a  6-degree  sector  an  x 
dimension  of  approximately  2  in.  would  be  produced  at  the  5-ft  range. 

However,  this  would  greatly  increase  the  number  of  elements  needed  for  the  SSI 
calculations.  Therefore,  the  effect  of  using  elements  with  variable  aspect 
ratios,  not  necessarily  close  to  one,  was  investigated.  Calculations  were 
performed  with  elements  using  a  constant  spacing  in  the  x  direction  and 
compared  to  those  with  a  constant  aspect  ratio  of  approximately  1.0. 
Calculations  were  also  performed  to  determine  the  optimum  constant  spacing. 

The  charge  in  the  tests  was  a  cylindrical  charge,  encased  in  steel,  containing 
15.4  lb  of  C-4.  A  sphere  containing  the  same  weight  of  C-4  has  a  radius  of  4 
in.  Therefore,  analyses  were  performed  using  a  4-in. -radius  sphere. 

40.  In  those  analyses  where  the  time  step  was  too  large  or  the  maximum 
strain  increment  or  viscosity  coefficients  were  too  low,  the  solution  was 
unstable  and  the  stress  time-histories  near  the  charge  look  similar  to  the  one 
shown  in  Figure  6.  Stable  solutions  could  not  be  obtained  using  the  default 
values  for  the  artificial  viscosity  coefficients.  In  some  cases  the  stress 
time-histories  further  away  appeared  to  be  acceptable.  Even  if  the  stresses 
further  away  from  the  charge  appear  to  be  acceptable,  the  parameters 
associated  with  these  analyses  were  rejected  because  the  effects  of  the 
instability  near  the  charge  cannot  be  evaluated  and  run  times  increased 
drastically  due  to  the  instability.  Once  a  stid>le  solution  was  obtained, 
further  parametric  calculations  were  performed  to  assess  the  accxiracy  the 
solution. 
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41 .  Figure  7  shows  a  comparison  of  the  results  of  computations  using  a 
constant  aspect  ratio  of  approximately  1  to  the  results  of  a  computation  using 
a  constant  element  thickness  of  1  in.  A  complete  description  of  these 
computations  is  provided  in  Table  2.  These  are  the  6-  and  2-degree  constant 
aspect  ratio  computations  described  above.  In  the  constant  aspect  ratio 
computation,  the  elements  are  thinner  near  the  charge  than  in  the  constant 
thickness  analyses.  This  results  in  higher  stresses  in  the  elements  near  the 
charge  and  higher  viscosities,  and  smaller  strain  increments  are  required. 
Figure  7  shows  that  the  rise  times  are  much  faster  in  the  constant  thickness 
computation.  This  is  because  near  the  points  of  interest,  the  elements  are 
much  larger  in  the  constant  aspect  ratio  computation,  and  the  rise  time  is 
clearly  a  function  of  the  element  thickness.  The  rise  time  also  increases  in 
the  constant  aspect  ratio  analyses  using  the  variable  radial  spacing  because 
of  the  increased  artificial  viscosity  required  in  these  runs.  Rise  times 
could  be  improved  by  taking  a  larger  number  of  elements  in  the  constant  aspect 
ratio  analysis,  but  this  would  make  the  elements  near  the  charge  extremely 
small,  and  since  the  critical  time-step  is  based  on  wave  travel  time  across  an 
element,  this  would  increase  computation  time  considerably.  The  rise  times 
associated  with  the  constant  thickness  elements  more  nearly  match  those  of  the 
test  data,  and  the  constant  thickness  element  was  selected  for  the  remaining 
calculations . 

42.  The  results  of  one  of  the  analyses  are  compared  with  experimental 
data  in  Figures  8  and  9.  Gage  locations  were  3,  4,  5,  6,  and  7  ft  from  the 
charge  with  the  stresses  (velocities)  arriving  at  the  gages  in  order,  based  on 
distance  from  the  charge.  In  this  analysis,  the  modified  Cap  parameters  were 
used.  The  viscosities  were  five  times  the  default  viscosities.  A  maximum 
strain  increment  of  0.0001  and  a  time-step  of  0.1  times  the  computed  critical 
time-step  were  used.  The  grid  consisted  of  1-in.  elements  up  to  40  ft  away 
from  the  center  of  the  charge.  These  figures  show  that,  qualitatively,  the 
results  compare  very  well  with  the  data.  However,  the  computed  peak  stresses 
and  velocities  were  much  lower  than  the  measured  values. 

43.  Since,  qualitatively,  the  results  were  good,  it  appears  that  the 
type  of  constitutive  model  being  used  was  acceptable.  A  number  of  parametric 
calculations  were  then  performed  in  an  attempt  to  tune  the  material  properties 
so  that  a  better  match  of  peak  stresses  and  velocities  with  the  data  could  be 
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obtained.  The  failure  surface  of  this  material  is  very  low,  and  essentially 
all  of  the  response  will  be  in  the  plastic  range.  Unless  the  failure  surface 
is  modified  drastically,  the  change  will  be  insignificant.  Therefore,  changes 
to  the  failure  surface  were  not  made.  The  bulk  modulus  and  hardening  function 
parameters  were  varied  to  determine  their  effect  on  the  stress  and  velocity 
time-histories.  These  runs  showed  that  changes  in  the  material  properties 
could  raise  the  stress  levels  to  close  to  rhe  measured  values  (Figure  10),  but 
the  peak  velocities  (Figure  11)  changed  very  little.  Figures  10  and  11  are 
the  stress  and  velocity  time-histories,  respectively,  at  5  ft  from  the  charge. 
In  this  analysis,  the  same  parameters  used  in  the  previous  analysis  were  used, 
except  for  the  maximum  volumetric  plastic  strain,  W,  of  0.003,  and  the 
hardening  function  parameter,  D1 ,  of  0.00138. 

44.  These  analyses  indicated  that  changes  in  material  properties  had 
very  little  effect  on  the  peak  velocities.  Since  the  characters  of  the  stress 
and  velocity  time-histories  were  a  good  match  to  the  data,  it  appears  that  the 
stress  wave  was  propagating  through  the  soil  correctly,  but  that  more  energy 
was  needed  from  the  explosive  source.  It  appears  that  the  JWL  model  is  not 
functioning  correctly.  This  will  be  discussed  more  fully  later  in  this 
report . 

45.  The  purpose  of  this  study  is  to  examine  the  effects  of  soil 
parameters  on  the  SSI,  not  to  predict  stresses  and  velocities  in  the  free 
field.  The  amount  of  energy  contributed  by  the  charge  can  be  increased  by 
simply  increasing  the  amount  of  charge  in  the  analysis.  Therefore,  analyses 
were  performed  using  the  modified  Cap  {>arameters  for  various  charge  sizes. 
Calculations  were  performed  for  5-,  6-,  7-,  8-,  and  9-in.  radius  charges.  It 
appears  that  a  charge  radius  between  6  and  7  in.  is  needed.  Therefore, 
analyses  were  performed  using  6-  and  7-in.  radius  charges  to  determine  the 
maximum  strain  increment,  critical  time-step  factor,  and  optimum  grid  sp>acing. 

46.  Preliminary  calculations  using  the  constitutive  driver  indicated 
that  at  a  maximum  strain  increment  above  0.0001,  the  model  becomes  unstable. 
Therefore,  this  was  the  highest  maximum  strain  increment  considered.  The 
DYNA3D  user's  manual  recommends  that  the  critical  time-step  factor  be  taken  as 
0.67  when  explosives  are  being  modeled.  Therefore,  this  is  the  largest  time- 
step  scale  factor  considered.  Calculations  were  performed  at  grid  spacings 
of  0.5,  1.0,  and  2.0  in. 
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47.  Figures  12  and  13  compare  the  stress  and  velocity  tire-histories, 
respectively,  at  a  range  of  5  ft  from  the  charge  for  calculations  using  1/2- 
and  1-in.  grid  spacings.  Experimental  data  are  also  included  for  comparison. 
These  figisres  show  that  the  rise  times  at  this  range  are  modeled  well  by 
either  the  1-in.  or  the  1/2-in.  grid.  It  appears  that  the  1-in.  grid  spacing 
will  be  adequate  since  there  is  very  little  difference  in  the  results  of  these 
two  runs.  A  comparison  between  the  results  using  a  2-in.  grid  spacing  and 
those  using  a  1/2-in.  spacing  is  shown  in  Figure  14.  This  figure  shows  that 
the  rise  tire  is  significantly  different  when  the  2-in.  grid  spacing  is  used. 
Therefoie,  a  1-in.  grid  spacing  was  selected  to  be  used  in  further 
calculations . 

48.  In  order  to  determine  the  importance  of  the  critical  time  step 
factor,  two  computations  were  performed.  Time-steps  of  0.1  and  0.67  times  the 
critical  time-step  were  used.  Radial  stress  and  velocity  time-histories  at 
the  five  different  ranges  were  compared,  and  the  results  of  the  two 
computations  were  identical;  therefore,  it  was  concluded  that  a  critical  time- 
step  factor  of  0.67  was  satisfactory. 

49.  Parametric  calculations  were  performed  for  maximum  strain 
increments  of  0.00001,  0.00002,  0.00005,  and  0.0001.  The  computations  using 
the  0.0001  increment  were  unstable.  The  radial  stress  time-histories  at  the 
3-'ft  range  for  the  other  three  analyses  are  compared  in  Figure  15.  This 
figure  shows  that  the  results  were  identical  for  the  0.00001  and  0.00002 
strain  increments,  while  the  arrival  time  was  slightly  greater  and  the  peak 
stress  was  less  for  the  analysis  using  the  0.00005  strain  increment.  The 
decays  of  the  stress  time-histories  were  approximately  the  same  for  all  three 
analyses.  Figure  16  shows  the  same  comparison  for  a  range  of  5  ft  from  the 
charge.  This  figure  also  shows  exact  agreement  between  the  computations  using 
strain  increments  of  0.00001  and  0.00002.  The  agreement  of  the  computations 
using  the  0.00005  strain  increment  with  the  other  analyses  was  much  better  at 
this  range.  The  arrival  was  only  slightly  later,  and  the  rise  time  and  peak 
stresses  were  very  nearly  the  same.  In  the  SSI  analyses,  the  structure  will 
be  placed  at  5  ft  from  the  charge.  Since  the  stress  time-history  predicted 
using  a  maximiun  strain  increment  of  0.00005  is  very  nearly  identical  to  those 
predicted  using  smaller  strain  increments,  this  strain  increment  should  be 
acceptable  for  performing  the  SSI  study. 
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50.  The  results  of  one  FE  analysis  are  compared  to  test  data  in  Figures 
17  and  18.  In  this  analysis  a  6-in.  charge,  1-in.  grid  spacing,  5  times  the 
default  viscosities,  a  critical  time-step  factor  of  0.67,  and  a  maximum  strain 
increment  of  0.00005  were  used.  Figxire  17  shows  that  analysis  results  compare 
extremely  well  with  the  stress  time-histories  at  the  3-  axid  4-ft  ranges,  but 
stresses  were  over  predicted  at  the  other  ranges.  At  the  5-,  6-,  and  7-ft 
ranges,  the  computed  rise  times  were  slightly  greater  and  the  stresses  decay 
slower  than  those  measured  in  the  test.  The  arrival  times  of  the  peak 
stresses  were  predicted  reasonably  well  at  all  ranges. 

51.  Figure  18  shows  that  the  velocity  time-histories  at  the  3-  and  4-£t 
ranges  were  predicted  well.  At  the  5-,  6-,  and  7-ft  ranges  the  velocities  are 
overpredicted.  Rise  times  were  predicted  reasonably  well  at  all  ranges,  but 
arrival  times  are  different  from  the  data  at  each  of  the  ranges.  The  arrival 
time,  based  on  velocity  at  a  given  range,  should  be  the  same  as  the  arrival 
time  based  on  stress  for  that  range.  This  is  true  based  on  the  analysis,  but 
is  not  true  based  on  the  data.  This  is  probably  due  to  the  methods  used  for 
collecting  the  velocity  and  stress  data. 

52.  These  analyses  required  that  relatively  high  artificial  viscosities 
be  used.  These  high  viscosities  cause  the  rise  times  of  the  stress  and 
velocity  time-histories  to  increase  with  distance  from  the  charge,  and  these 
rise  times  were  too  long  as  compared  to  test  data.  Axialyses  were  performed  to 
determine  if  using  artificial  viscosity  coefficients  that  varied  with  distance 
from  the  charge  would  improve  the  rise  times  at  the  ranges  farther  away  from 
the  charge.  An  analysis  identical  to  the  one  shown  previously  was  performed, 
except  that  the  maximum  strain  increment  was  0.00002,  and  variable  artificial 
viscosity  terns  were  used.  When  variable  artificial  viscosities  were  used, 
there  was  a  significant  difference  between  the  computation  using  a  maximum 
strain  increment  of  0.00002  and  the  one  using  0.00005.  Therefore,  0.00002  was 
selected.  In  this  analysis,  the  following  multiples  of  the  default 
viscosities  were  used: 
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Range  from  charge,  in.  Multiple  of  default  viscosities 

5 
4 
3 
2 
1 

53.  Figure  19  compares  the  results  of  this  analysis  with  the  previous 
analysis  and  test  data.  There  were  no  differences  between  the  stress  or 
velocity  time-histories  for  these  two  analyses  at  the  3-  and  4-£t  ranges; 
therefore,  those  ranges  are  not  shown.  At  the  5-ft  range,  there  was  a  slight 
improvement  in  the  rise  time.  At  the  6-  and  7-ft  ranges,  there  was  a 
significant  improvement  in  both  rise  times  and  arrival  times.  Figure  20  shows 
that  the  arrival  times  of  the  velocities  did  not  agree  as  well  with  the 
velocity  data,  but  the  rise  times  using  the  variable  viscosities  agreed  better 
with  the  data  at  the  5-,  6-,  and  7-ft  ranges.  These  analyses  showed  that 
using  artificial  viscosities,  which  vary  with  range,  can  significantly  improve 
the  rise  times  of  the  stress  and  velocity  time-histories.  The  use  of  these 
variable  viscosities  will  produce  stress  and  velocity  time-histories  which  are 
more  like  those  measured  in  the  test. 

Cylindrical  Clav  Analyses 

54.  Since  a  2-D  grid  may  be  used  for  many  of  the  SSI  studies,  the 
effects  of  using  a  cylindrical  geometry  as  opposed  to  a  spherical  geometry 
should  be  investigated.  It  is  desirable  that  the  correct  stress  and  velocity 
time-histories  be  predicted  at  least  at  the  structure  location.  Therefore, 
the  objective  of  these  analyses  is  to  determine  if  a  charge  size  can  be 
selected  so  that  the  correct  stress  and  velocity  time-histories  are  predicted 
at  the  structure  location.  Since  the  stress  and  velocity  gradients  at  the 
structiire  will  probably  affect  the  SSI  problem,  the  effect  of  the  2-D  geometry 
on  stress  gradients  should  also  be  Investigated,  ideally,  near  the  structure 


IR  <  42 

42  <  R  <  48 
48  <  R  <  60 
60  <  R  <  72 
72  <  R 
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the  attenuation  of  peak  stress  (velocity)  versus  range  should  be  the  same  in 
order  for  the  stress  (velocity)  gradients  to  be  approximately  the  same. 

55.  VE  analyses  were  performed  using  a  cylindrical  charge  with  a  1>in. 
radius.  The  computations  were  performed  using  a  1-D  grid  similar  to  the  one 
used  for  the  spherical  analyses.  In  this  case,  the  elements  were  1-in.  thick 
in  the  z  direction  and  grew  in  the  y  direction  as  x  increased.  A  6-degree 
sector  of  a  cylinder  was  used  for  these  calculations,  and  boundary  conditions 
were  used  to  enforce  cylindrical  symmetry.  Material  properties  were  the  same 
as  those  used  for  the  spherical  calculations.  In  these  analyses,  a  critical 
time-step  factor  of  0.25  and  a  maximum  strain  increment  of  0.00002  were  used. 
The  artificial  viscosity  factors  were  as  followi;: 

Range  from  charge,  in.  Multiple  of  default  viscosities 

R  <  6  20 

6  <  R  <  12  15 

12  <  R  <  24  10 

24  <  R  <  36  8 

36  <  R  <  48  5 

48  <  R  3 


56.  Figure  21  shows  that  the  stresses  at  the  5-ft  location  in  the 
cylindrical  analysis  agree  very  well  with  those  of  the  spherical  analysis. 
Figure  22  shows  that  the  maximum  velocities  agree  reasonably  well  between  the 
two  analyses,  but  the  velocities  decay  slightly  faster  for  the  cylindrical 
analyses.  Peak  stresses  versus  range  from  the  charge  are  shown  in  Figure  23. 
This  figure  shows  that  near  the  5-ft  range  the  stresses  in  the  cylindrical 
analyses  attenuate  with  range  similarly  to  those  in  the  spherical  analysis, 
and  both  analyses  match  the  data  reasonably  well.  Figure  24.  shows  the 
attenuation  of  maximum  velocity  with  range.  Tliis  figure  shows  that  the 
cylindrical  analysis  predicts  the  attenuation  of  maximum  velocity  with  range 
reasonably  well.  These  figures  indicate  that  it  is  reasonable  to  use  a  2-D 
model  to  perform  the  SSI  studies. 

57.  Since  artificial  boundaries  must  be  introduced  into  the  SSI 
analysis,  it  is  important  to  determine  how  well  the  nonreflecting  boundaries 
work.  In  preliminary  studies  to  eval\iate  the  noiureflecting  boundaries,  it  was 
discovered  that  the  nonreflecting  boundary  was  not  affected  by  changing  the 
properties  of  the  material  near  the  boundary.  In  investigating  this  problem 
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it  wa8  determinod  that  the  subroutine  which  passes  material  property 
information  to  the  nonreflecting  boundary  subroutine  was  not  passing  the 
correct  information.  This  vjas  corrected  and  further  studies  were  performed. 
Figures  25  and  26  show  the  effects  of  different  boundary  conditions  on 
stresses  and  velocities,  respectively,  at  5  ft  from  the  charge.  The  solid 
lines  in  these  figures  are  for  a  nonreflecting  boundary  at  40  ft.  This 
boundary  was  far  enough  away  so  that  boundary  effects  did  not  appear  for  the 
time  shown  in  these  figures.  The  other  results  are  for  various  types  of 
boundaries  placed  at  10  ft  from  the  charge.  These  figures  show  that  the 
stresses  and  velocities  based  on  the  nonreflecting  boundary  at  10  ft  were  very 
close  to  those  using  the  fixed  boundary.  This  is  not  satisfactory  for  the  SSI 
study . 

58.  The  nonreflecting  boundary  is  designed  to  simulate  an  infinite 
amount  of  material  beyond  the  boundary.  In  DYNa3D  the  model  is  based  on 
constant  values  of  the  shear  and  bulk  moduli,  and  is,  thus,  only  accurate  if 
the  boundary  is  placed  in  a  region  such  that  the  soil  is  elastic.  This  is  not 
true  in  this  case,  and  the  error  introduced  is  obviously  significant.  Figure 
27  shows  a  comparison  of  several  forms  of  the  volumetric  stress-strain  curve 
for  this  soil.  Clearly,  using  the  elastic  constants  only  will  produce  a 
boundary  that  is  much  stiffer  than  it  should  be.  The  material  model  was 
modified  so  that  the  shear  and  bulk  moduli  used  by  the  nonreflecting  bouxidary 
subroutine  could  be  different  from  those  used  by  the  constitutive  model  for 
those  elements. 

59.  The  volumetric  stress-strain  curve  needed  to  tune  the  nonreflecting 
boundary  is  shown  in  Figure  27 .  This  shows  that  the  elastic  consta^its  must  be 
selected  so  that  the  elastic  stress-strain  curve  for  the  boundary  matches  the 
total  stress-strain  curve  for  the  material.  The  elastic  bulk  and  shear  moduli 
were  each  divided  by  80  in  this  case. 

60.  Figures  28  and  29  show  comparisons  of  stress  and  velocity 
time-histories,  respectively,  based  on  analyses  using  the  modified 
nonreflecting  boundary  at  10  ft  from  the  charge,  to  those  in  which  the 
boundary  is  too  far  away  to  be  significant.  These  figxjres  show  that  the 
modified  boundary  simulates  a  nonreflecting  boundary  very  well. 
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61 .  The  material  properties  of  sand  are  significantly  different  from 
those  of  clay.  Figure  30  shows  that  the  failure  surface  of  clay  is  almost 
constant  and  is  much  lower  than  that  of  aax^.  The  failxire  surface  in  sand 
increases  significantly  with  increasing  mean  normal  stress.  Figure  31  shows  a 
comparison  of  the  volumetric  stress-strain  curves  for  the  clay  and  sand  tested 
in  the  CONWEB  test  series.  This  shows  that  the  stress-strain  curve  of  the 
clay  locks  up  at  approximately  4  percent,  while  that  of  the  sand  locks  up  at 
approximately  26  percent.  These  numbers  correspond  to  the  percentage  of  the 
volume  of  each  material  which  is  initially  filled  with  air,  and  indicate  that 
once  these  air  voids  are  closed  further  plastic  volumetric  straining  does  not 
occur. 

62.  One-dimensional  spherical  analyses  similar  to  those  performed  for 
the  clay  material  were  performed  for  the  saxui  material.  Material  properties 
for  the  Cap  model  for  the  sand  were  based  on  static  tests.  These  material 
properties  are  listed  in  Table  1 .  As  was  the  case  with  the  clay,  parametric 
calculations  were  performed  to  determine  the  required  grid  spacing,  critical 
time-step  ratio,  artificial  viscosity  terms,  and  naximxim  strain  increment. 

63.  In  crder  to  obtain  suitable  results,  the  value  of  the  shape 
parameter  D1  of  the  hardening  function  was  changed  from  0.0000758  per  psi  to 
0.00003  per  psi.  All  other  material  properties  were  as  listed  in  Table  1. 

The  material  properties  using  this  value  of  D1  will  be  referred  to  as  the 
modified  properties  for  sand.  The  effect  of  this  change  in  D1  on  the 
volumetric  stress-strain  cxirve  is  shown  in  Figure  31 .  The  following 
parameters  were  needed  to  reasonably  predict  the  stresses  and  strains  in  the 
sand  test: 

Charge  Radius  6  in. 

Grid  Spacing  1  in. 

Time-Step  Ratio  0.67 

Artificial  Viscosity  Terms  3  times  default 

Maximum  Strain  Increment  0.00001 

64.  Figures  32  and  33  compare  the  stresses  and  velocities, 
respectively,  from  this  analysis  with  the  test  data.  The  computed  arrival 
time  of  the  stresses  at  3  ft  did  not  agree  with  the  data;  therefore,  0.3  msec 


was  subtracted  from  all  of  the  data  records  to  account  for  this  difference. 
Figure  32  shows  that  stresses  were  predicted  reasonably  well  for  this  test. 

In  general,  the  arrival  times  were  sooner  in  the  calculation,  and  the  rise 
times  were  longer.  The  arrival  times  become  progressively  worse  with  distance 
from  the  charge.  For  those  measurements  at  3,  4  and  5  ft,  the  analysis 
predicted  that  the  stresses  will  drop  more  quickly  from  their  peak  values  than 
was  measured  in  the  test.  Late-time  stresses  were  significantly  overpredicted 
at  all  ranges  as  shown. 

65.  Figure  33  shows  that  the  comparison  of  measured  to  computed  free- 
field  velocities  is  similar  to  the  comparison  of  stresses.  In  general,  the 
predicted  maximum  velocities  agreed  well  with  the  test  data.  The  predicted 
arrival  times  were  too  soon,  and  the  rise  times  were  too  long,  unfortunately 
two  of  the  gages  (4  and  5  ft)  malfunctioned  after  peak  velocity  was  obtained, 
and  complete  velocity  time-histories  were  not  available.  The  data  for  the  7- 
ft  range  are  not  consistent  with  the  remainder  of  the  velocity  data  for  this 
test,  and  this  comparison  with  the  test  data  cannot  be  made.  For  the  other 
two  gages,  late-time  free-field  velocities  were  underpredicted  .  This  is 
consistent  with  the  overprediction  of  late-time  free-field  stresses  and 
indicates  that  the  free  surface,  which  was  not  considered  in  the  analysis, 
affected  both  late- time  stresses  and  velocities  in  the  experiment  in  sand. 

The  relatively  small  and  shallow  test  bed  of  sand  was  also  surrounded  by  the 
native  clay  backfill  material,  similar  to  the  clay  previously  analyzed.  This 
effect,  which  was  not  modeled  in  the  auialysis,  would  also  affect  late-time 
stresses  and  velocities. 

66.  Fiirther  analyses  were  performed  using  artificial  viscosity  terms 
which  varied  with  distance  from  the  charge  to  determine  if  rise  times  and 
arrival  times  would  be  improved.  The  following  multiples  of  the  artificial 
viscosity  terms  were  used: 

Range  from  charge  center,  in.  Multiple  of  default  viscosities 
R  <  9  3 

9  <  R  <  13  2 

13  <  R  <  18  1 

18  <  R  <  24  0.5 

24  <  R  0.25 
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67.  Figure  34  shows  that  arrival  times  and  rise  times  of  free-field 
stresses  were  significantly  improved  at  all  ranges  when  variable  artificial 
viscosities  were  used.  The  histories  were  shifted  so  that  the  peak  stresses 
at  3  ft  would  occur  at  the  same  time  for  the  two  analyses  and  the  data.  The 
peak  stresses  were  increased  slightly  in  the  analysis  using  the  variable 
artificial  viscosities;  and  since  the  peak  stresses  were  already  slightly 
high,  the  agreement  between  predicted  and  measured  peak  stresses  was  worsened. 
Agreement  with  late-time  stresses  was  not  changed  at  all.  Since  the  arrival 
and  rise  times  were  significantly  improved  with  only  a  small  increase  in  error 
in  peak  stresses,  the  overall  effect  of  using  the  variable  artificial 
viscosities  was  to  Improve  the  agreement  between  the  predicted  and  measured 
stress  time-histories. 

68.  The  effect  of  variable  artificial  viscosity  on  free-field 
velocities  is  shown  in  Figure  35.  Arrival  times  and  rise  times  were  both 
significantly  improved  by  using  the  variable  artificial  viscosity.  Maximum 
predicted  free-field  velocities  were  higher  for  the  variable  viscosity 
calculation,  and  agreement  with  the  data  was  slightly  worsened.  The  overall 
effect  of  using  the  variable  artificial  viscosities  was  to  improve  the 
agreement  of  the  computed  velocity  time-histories  with  the  test  data. 

Cylindrical  Sand  Analyses 

69.  Analyses  were  performed  to  determine  the  effects  of  using 
cylindrical  versus  spherical  boundary  conditions  in  the  sand  backfill 
material.  A  1-in.  radius  charge  was  used  in  this  study.  Parametric 
calculations  were  performed  until  a  good  combination  of  the  critical  time 
step,  maximum  strain  increment,  and  artificial  viscosity  terms  was  determined. 
The  analyses  were  performed  using  the  same  grid  as  that  used  for  the 
cylindrical  analyses  in  clay.  The  material  properties  were  the  same  as  those 
used  for  the  spherical  sand  analyses.  In  the  cylindrical  sand  analyses,  a 
critical  time-step  factor  of  0.25  and  a  oiaximum  strain  increment  of  0.00001 
were  needed.  The  following  multiples  of  the  default  artificial  viscosity 
terms  were  used: 
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Multiple  of  default  viscosity 
30 


Radivis  from  center  of  charge,  in. 


R  <  3 

3  <  R  <  5  20 

5  <  R  <  9  10 

9  <  R  <  16  5 

16  <  R  <  24  3 

24  <  R  <  46  2 

46  <  R  1 


70.  Figures  36  and  37  show  comparisons  between  the  analyses  using 
cylindrical  versus  spherical  boundary  conditions  for  stress  and  velocity  time- 
histories,  respectively,  for  a  range  of  5  ft  from  the  center  of  the  charge. 
These  figures  show  that  the  1-in.  charge  did  a  good  job  of  matching  the  stress 
and  velocity  time-histories  at  a  range  of  5  ft.  Both  stresses  and  velocities 
were  slightly  low.  Indicating  that  there  probably  is  a  slightly  larger  charge 
that  would  do  a  better  job  of  matching  the  peak  stresses  and  velocities. 
However,  this  analysis  does  indicate  that  an  analysis  imposing  cylindrical 
boundary  condition  can  be  used  to  match  the  stress  and  velocity  time-histories 
at  the  5-ft  range  from  a  spherical  source. 

71 .  Figures  33  aivd  39  show  the  attenuations  versus  range  of  peak  stress 
and  maximum  velocity,  respectively,  for  the  sidierical  and  cylindrical 
geometries.  These  figtures  show  that  the  attenuations  versus  range  of  peak 
stress  and  maximxim  velocity  of  the  cylindrical  analysis  match  those  of  the 
spherical  analysis  reasonably  well  near  the  5-ft  range.  Either  analysis 
matches  the  data  reasonably  well.  Based  on  Figures  36-39  it  is  reasonable  to 
use  a  2-D  model  to  study  SSI  in  the  sand  material. 

72.  Figures  40  and  41  show  the  comparison  of  stresses  and  velocities, 
respectively,  at  the  5-ft  range  from  the  charge,  between  an  analysis  with  the 
boundary  far  away,  and  another  analysis  in  which  a  modified  nonreflecting 
boundary  was  placed  at  10  ft  from  the  charge.  The  boundary  was  modified  by 
dividing  the  elastic  bulk  and  shear  moduli,  respectively,  by  50.  Figures  40 
and  41  show  that  the  modified  nonreflecting  boundary  did  an  excellcmt  job  for 
at  least  20  msec. 
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73.  It  was  noted,  in  a  previous  section  that  the  charge  size  needed  in 
the  spherical  calculations  to  match  the  data  from  the  CONWEB  experiments  was  a 
6-in. -radius  charge,  while  the  charge  used  in  the  experiment  was  equivalent  to 
a  4-in. -radius  charge.  Although  it  was  not  necessary  for  this  study,  it  was 
desirable  to  understand  the  reason  for  this  discrepancy. 

74.  In  these  analyses,  two  tetrahedral  elements  were  used  to  model  the 
explosive.  In  an  effort  to  determine  the  source  of  this  error,  it  was 
discovered  that  the  volume  of  the  tetrahedrons  used  to  model  the  explosives 
was  computed  Incorrectly.  The  soxirce  of  this  error  is  discussed  in  Appendix 
A.  This  error  does  not  significantly  affect  the  JWL  model,  but  does 
significantly  affect  the  nodal  loads  computed  based  on  the  pressure  in  the 
element.  For  the  6-in.-radxus  charge  a  volume  of  0.29  cu  in.  was  computed, 
while  the  actual  volume  was  approximately  0.78  cu  in. 

75.  If  a  very  small  portion  of  the  volume  is  deleted  from  the  center  of 
the  sphere,  the  sphere  can  be  modeled  using  an  8-node  solid  element.  Although 
this  element  has  an  extremely  bad  aspect  ratio  (approximately  20  based  on  the 
small  dimension  near  the  center  of  the  element),  the  computed  volume  is  very 
nearly  correct.  Analyses  were  performed  for  the  sand  material  using  a 

4-in. -radius  charge  with  0.2  in.  removed  from  the  center.  Thus,  the  cliarge 
can  be  isodeled  using  a  single  8-node  element,  and  the  volume  is  computed  very 
nearly  correctly  (14  percent  low).  In  these  analyses,  the  modified  sand 
properties  were  used.  The  critical  time-step  factor,  artificial  viscosity 
factor,  and  maximum  strain  increment  were  0.67,  20,  and  0.00001,  respectively. 
The  artificial  viscosity  factor  was  dropp>ed  to  5  at  a  distance  of  1^  in.  from 
the  center  of  the  charge. 

76.  Figures  42  and  43  compare  the  stress  and  velocity  time-histories, 
respectively,  from  this  analysis  to  the  test  data.  These  figures  show  that 
the  stress  and  velocity  time-histories  compare  reasonably  well  with  the  test 
data.  This  indicates  that  a  6- in. -radius  charge  was  needed  in  the  analysis 
using  the  tetrahedral  elements  for  the  charge  because  the  volume  of  the  charge 
was  computed  incorrectly,  and  that  more  care  should  be  talien  in  developing  the 
grid  for  the  explosives.  Thi  volume  of  the  6-node  element  used  in  the 
cylindrical  calculations  was  computed  correctly  by  DYNA3D. 
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PART  XV:  SUMMARY  AND  CONCLUSIONS 

SwRffiary 

77.  FE  analyses  which  lncl^lde  the  detonation  of  the  explosive  charge, 
propagation  of  stresses  through  the  soil  to  the  structure,  interaction  of  the 
soil  with  the  structure,  and  response  of  the  structure  are  needed  to  study 
SSI.  DYNA3D  can  model  the  detonation  of  the  explosive  using  the  JWL  equation 
of  state.  It  also  contains  a  version  of  the  Cap  model,  which  has  been 
demonstrated  to  be  very  effective  for  modeling  wave  propagation  through  soil 
materials.  Subroutines  in  the  code  are  available  for  modeling  the  interaction 
of  the  soil  with  the  structvure,  and  material  and  geometric  nonlinearities  can 
be  analyzed.  Therefore,  0YNA3D  appears  to  be  ideal  for  investigative  SSI. 
OYNASD  is  widely  available  to  engineers  who  would  like  to  perform  these 
calculations,  and  prior  to  this  stiidy  there  were  no  published  reports  of  SSI 
analyses  of  this  typo  using  DYNA3D  or  any  other  FE  code.  It  was  therefore 
assumed  that  using  DYNA3D  for  the  SSI  investigation  will  be  difficult,  if  not 
impossible.  This  study  was  performed  to  determine  if  the  detonation  of  the 
explosive  and  propagation  of  the  stress  wave  through  the  soil  could  be 
adequately  modeled. 

78.  Preliminary  analyses  showed  that  the  Cap  model  in  DYNASD  would  not 
function  for  the  material  properties  to  be  uscsd  in  this  study.  Another  Cap 
model  was  obtained  and  installed  in  DYNA3D.  This  Cap  model  uses  a  constant 
bulk  modulus.  This  is  adequate  for  the  soil  types  being  analyzed  as  long  as 
pressures  are  not  extremely  high.  However,  at  locations  near  the  charge  a 
bulk  modulus  which  increases  with  increasing  pressure  is  needed.  The  newly 
installed  Cap  model  was  modified  to  use  a  bulk  modulus  which  increases  with 
increasing  pressure.  This  Cap  model  was  also  modified  so  that  the 

strain  increment  is  input  instead  of  being  computed  by  the  code.  This 
the  Cap  model  much  more  efficient. 

79.  One-dimensional  spherical  calculations  for  comparison  to  the  CONffEB 
clay  and  sand  tests  were  then  attempted.  In  these  calculations,  the  DYNA30 
default  values  of  the  viscosity  coefficients  were  used  and  the  runs  were 
unstable.  It  was  determined  that  much  higher  viscosity  coefficients  were 
needed  to  stabilize  the  results.  However,  when  these  large  viscosities  were 
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used,  the  rise  times  at  locations  away  from  the  charge  were  much  too  high. 

When  viscosity  coefficients  which  decreased  with  increasing  distance  from  the 
charge  were  used,  the  calculations  were  stable  and  the  rise  times  were 
acceptable . 

80.  In  adoition,  optimum  values  of  grid  spacing,  maximum  strain 
increment,  and  critical  time  step  for  the  sand  and  clay  materials  were 
investigated.  The  results  of  there  analyses  compared  reasonably  well  with 
experimental  data. 

81 .  In  order  to  determine  if  the  SSI  study  could  be  conducted  in  2-D, 
one-dimensional  cylindrical  calculations  were  performed  and  compared  to  the 
1-D  spherical  calculations.  At  the  range  of  the  structure,  the  same  peak 
stresses  and  velocities  as  predicted  for  the  spherical  calculations  could  be 
obtained  by  using  a  cylindricil  geometry  and  a  smaller  charge  radius.  The 
shapes  of  the  stress  and  locity  time-histories  at  5  ft  were  also  reasonably 
close  to  those  of  the  s- her leal  calculations.  At  the  structure  range,  peak 
stress  and  maximum  velocity  gradients  in  the  cylindrical  calculation  also 
agreed  very  well  with  the  spherical  calculation. 

82.  The  properties  needed  for  the  nonreflecting  boundary  subroutine 
were  not  being  passed  to  that  subroutine.  This  was  corrected,  however,  the 
noiureflecting  boundary  was  still  much  too  stiff  blnce  the  routine  is  based  on 
the  elastic  bulk  and  shear  moduli  and  is,  therefore,  too  stiff  if  a 
nonreflecting  boundary  is  placed  in  an  inelastic  region.  The  subroutine  was 
modified  so  that  a  fraction  of  the  elastic  properties  could  be  used.  It  was 
demonstrated  that  if  the  correct  fraction  of  the  elastic  material  properties 
is  used,  the  nonreflecting  boundary  method  in  DYNASD  does  a  good  job  of 
simulating  the  continuum  beyond  the  boundary.  Factors  for  a  boundary  at  10  ft 
from  the  charge  were  determined  for  the  clay  and  sand  materials. 

Conclusions 

83.  The  Cap  model  in  DYNA3D  did  not  function  correctly  for  the  sand  and 
clay  backfill  materials  used  in  the  CQNWEB  experiments.  Stable  solutions 
could  not  be  obtained  using  the  default  viscosities  and  the  nonreflecting 
boundaries  in  DYNA3D.  Therefore,  the  vnaltered  version  of  DyNA3D  could  not  be 
used  to  investigate  SSI.  The  new  Cap  model  did  function  well  for  the  clay  and 
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sand  backfills  and  when  it  was  modified  to  use  a  variable  bulk  modulus  it 
could  be  used  to  model  the  soil  adjacent  to  the  charge.  The  modified 
nonreflecting  boundary  worked  extremely  well.  DYNA3D  using  the  modified  Cap 
model,  variable  artificial  viscosities,  and  the  modified  nonreflecting 
boundary  is  very  suitable  for  performing  the  SSI  study.  It  was  also 
determined  that  the  SSI  study  could  be  conducted  in  2-D. 
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Table  1 


Properties  of  Clav  and  Sand 


tt 

Parameter 

units 

Sssd 

psi 

1.16  X  10® 

2.219  X  10® 

psi 

312 

312 

K2 

0.7 

0.7 

G 

psi 

43,511 

1.088  X  10® 

a 

psi 

16.68 

18,564.69 

e 

0.0 

0.0 

B 

1/psi 

0.00434 

2.530  X  10’ 

Y 

psi 

13.78 

18,562.51 

R 

2.5 

3.5 

D1 

1/psi 

0.00138 

7.585  X  10" 

02 

1 /psi^ 

0.00 

0.00 

W 

0.04 

0.256 

psi 

0.00 

0.00 

p  Ib- 

2/  •  4 

■sec  /in. 

1.841  X  10”^ 

1.752  X  lO"" 

T 

psi 

0.00 

0.00 

Table  2 

Parameters  for  Grid  Spacing  Rims 


Run 

G-degroe  constant 
aspect  ratio 
2 -degree  constant 
aspect  ratio 
1-in.  constant 
radial  spacing 


*Max 

viscosity 

0.00001 

10 

10 

0.00001 

10 

10 

0.0001 

10 

5 

Parameters  are  defined  in  paragraphs  24-26. 
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Figure  1.  Failure  stirface  fcr  cap  aiodel  in  DYMA3I 


Figure  2.  Hardening  function  in  DYNABD  cap  model 


Figure  3.  Comparison  of  hydrostatic  stress-strain  curves  for  clay 


Figure  4.  One-dimensional  spherical  grid 
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Figure  5.  Four-node  element  for  charge 
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r«i  14.  Effect  of  1/2“in.  vc  2-in  grid  spacing  on  stress  time-history 
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Figure  15.  Effect  of  changes  in  maximum  strain  increment  on  stresses  at  3  ft 
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Figure  17.  Comparison  of  measured  to  computed  radial  stresses  (continued) 
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Figure  18.  Comparison  of  measured  to  computed  radial  velocities  (continued) 
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Figure  18.  (Concluded) 


^***‘  P*'  fitfMS.  psi  efr,,,. 


a.  Range  •■  5ft 


b.  Range  6ft 


Time,  msec 
c  Range  -  7  ft 


Figure  19.  l!££ect  of  v^iriable  artificial  viscosity  on  stresses 
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Figurtt  23.  Comparifion  of  cylindrical  to  apharical  (paak  atrasa  varaua  ranga) 


Figure  24.  Comparison  of  cylindrical  to  spherical  (peak  velocity  versus  range) 
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Figure  25.  Effect  of  boundary  conditions,  stresses  at  5  ft 


Figure  26.  Effect  of  boundary  conditions,  velocities  at  5  ft 
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Figure  31 .  Comparison  of  Hydrostatic  stress-strain  curves  for  clay  euid  sand 
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Figure  33.  Comp^u:ison  of  measured  to  computed  radial  velocities 

for  sand  (continued} 
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Figure  33.  (Concluded) 
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55  Effect  of  variable  viscosits  on  velocities 
Figure  35.  ^  (continued) 
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Figure  38.  Comparison  of  cylindrical  to  spherical  in  stand 

peak  stress  versus  mage 
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Figure  39.  Comparison  of  cylindrical  to  spherical  in  sand 
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Figure  40.  Effect  of  modified  nonreflecting  boundary 
on  Gtresses  at  5  ft  in  sand 


Figinre  41.  Effect  of  modified  nonreflecting  boundary 
on  velocity  at  5  ft  in  sand 
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Figiire  42.  Stresses  computed  using  8-node  element 
for  charge  (continued) 


■ 


RMfial8trMs,p«l  ftodtol  Strm.  ptf 


400 


TInw,  msac 
•.  Rang*  b  7  ft 


Figure  42.  (Concluded) 
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APPENDIX  A:  CALCULATION  OF  VOLUME  OF  SOLID  ELEMENTS  IN  DYNA3D 


1 .  In  DYNABD,  the  time-step  for  a  given  element  is  determined  based  on 
the  dimensions  of  that  element  and  the  properties  of  the  material  of  which 
that  element  is  composed.  The  critical  time-step  is  the  smallest  tine-step  of 
all  of  the  elements  making  up  the  grid.  If  the  charge  is  broken  up  into  a 
large  number  of  elements/  the  elements  and,  thus,  the  critical  time-step  would 
be  very  small.  This  small  time-step  would  control  the  time-step  for  the  FE 
analysis.  If  larger  elements  could  be  used  for  the  charge,  the  critical  time- 
step  would  increase  and  the  analyses  could  be  performed  more  efficiently. 

2.  The  explosive  is  modeled  using  the  JWL  equation  of  state,  which  is 
accurate  as  long  as  the  initial  and  current  volume  of  the  element  are  computed 
correctly.  Since  the  volume  can  '^ery  easily  be  computed  exactly,  even  for  a 
severely  distorted  element,  it  was  decided  rhat  the  explosive  could  be  modeled 
using  only  one  element  in  the  radial  direction.  For  the  1-D  spherical 
analyses,  the  charge  was  modeled  using  two  tetrahedral  elements. 

3.  The  solid  elements  used  in  DYNA3D  are  constant  strain  trilinear 
hexahedral  elements,  as  discussed  in  Hughes  (1987).  The  DYNA3D  user's  manual 
(Hailquist  and  Benson  1987)  states  that  4-,  6-,  and  8-node  elements,  such  as 
those  shown  in  Figure  Al  may  be  used.  This  is  an  isoparametric  element  where 
the  parent  element  (Figure  Al )  is  mapped  into  the  actual  element  geometry, 
which  may  have  between  four  and  eight  nodes.  The  parent  elesment  is  a  cube 

^ts  center  at  the  origin  of  the  coordinate  system,  A.  The  coordinate  of 
each  of  the  corners  ds  either  plus  or  minus  1.  For  example,  node  1  is  at  (-1, 
-1,  -1),  and  node  7  i.,  at  (1,  1,  1). 

4.  The  parent  element  is  mapped  into  the  solid  element  using  the 
following  relationships: 


x(a,  b,  c) 


y(a,  b,  c) 


z(a,  b,  c) 
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where  x^,  y^,  z.  are  the  x,  y,  and  z  coordinates  of  the  ith  node,  and  is 

Al 


the  ith  shape  function  given  by: 

H^(a,  b,  c)  =  1/8(1  +  a^a){1  +  b^b) (1  +  c^c) 

5 .  These  same  functions  are  also  used  to  compute  displacements  in  the 
parent  element  from  the  displacements  of  the  node  points,  for  example: 

u(a,  b.  c)  ■> 


where  u  is  the  displacement  in  the  x  direction 
u^  is  the  X  displacement  of  the  ith  node 
The  Jacobian  of  the  transformation  from  the  A  system  to  the  X  system  is  given 
by: 

J  =  y.,  y,t  y'c 

*'a  *»b  *'c 

k  J 

where  the  comma  indicates  a  partial  derivative  of  the  variable 
with  respect  to  the  subscript  variable.  For  example: 
is  the  partial  derivative  of  x  with  respect  to  a. 

This  Jacobian  matrix,  or  its  inverse,  is  used  to  compute  nodal  loads  and 
strains  in  the  element.  Since  this  is  a  constant  strain  element,  the  strains 
may  be  computed  at  the  origin  of  the  patent  element. 

6.  The  original  volume  of  the  element  is  computed  as  the  integral  over 
the  volume  of  the  parent  element  of  the  determiniuit  of  the  Jacobian  matrix. 

The  current  volume  may  be  computed  using  the  current  coordinates  of  the  nodes 
in  computing  the  Jacobian  matrix.  Consistent  with  using  tlie  strains  at  the 
origin  of  the  parent  system,  the  volume  ca^i  be  approximated  using  a  single 
Gauss  point  at  the  origin.  Thus,  the  volume  is  approximated  as  S  times  the 
value  of  the  determinant  of  the  Jacobian  matrix,  evaluated  at  the  origin.  The 
volume  is  computed  exactly  if  the  determinant  of  the  Jacobian  is  constant  or 
an  odd  function  of  the  parent  coordinates.  If  the  determinant  is  neither  a 
consteint  nor  an  odd  function  of  the  parent  coordinates,  the  volume  is  only 
approximated. 

7.  For  tne  elements  used  to  model  a  6-degree  sector  of  a  sphere  with  a 
radius  of  4.0  in.,  the  volme  computed  from  the  values  cf  the  determinants  at 
the  origins  of  the  two  elements  is  0.08720  in.  The  actual  total  volume  of  the 
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two  elements  is  0.23180.  Thus,  the  computed  volume  of  the  explosives  is 
approximately  60  percent  low. 

8.  Because  the  nodal  loads  are  computed  using  similar  procedures  as 
used  to  compute  the  volume,  the  nodal  loads  are  also  computed  incorrectly. 

This  causes  the  deformations  of  the  grid  to  deviate  from  spherical  symmetry. 
The  boundary  ccndj  tions  ensure  that  the  nodes  move  along  radial  lines  through 
the  center  of  the  charge  but  do  not  ensure  that  each  node  at  a  given  radius 
moves  the  same  amount.  Apparently  the  nodes  shared  by  the  elements  receive 
more  load  than  the  other  two,  since  those  shared  nodes  move  faster.  This 
difference  in  '’elocity  is  noticeable  only  in  the  nodes  near  to  the  charge,  and 
there  is  very  little  difference  at  distances  of  interest  in  this  study. 

9.  Given  the  method  of  computing  volumes  that  is  used  in  DYNA3D,  it  is 
important  to  be  careful  in  developing  the  grid  to  be  used  for  the  explosives 
in  these  calculations.  The  error  in  volume  can  be  reduced  in  several  ways.  In 
the  calculations  summarized  in  Part  3,  a  charge  with  a  6-in.  radius  was  used. 
The  additional  volume  caused  by  the  increase  in  radius  offsets  the  fact  that 
the  volume  is  undercomputed  in  the  code.  In  this  study  the  objective  was  to 
demonstrate  that  this  method  could  be  used  to  compute  stresses  and  velocities 
at  various  ranges.  Using  the  6-in.  charge  is  an  acceptable  way  to  demonstrate 
this. 

10.  Another  way  to  minimize  the  volume  error  is  to  use  a  5-node  solid 
for  the  explosive.  Even  though  the  user's  manual  indicates  that  this  element 
is  not  available,  there  is  no  reason  why  this  element  can't  be  used.  Although 
no  satisfactory  calculations  were  performed  using  this  element,  the  volume 
computed  using  DYNA3D  is  only  14  percent  lower  than  the  actual  volume  of  the 
explosives,  and  displacements  are  spherically  symmetric. 

11.  A  third  way  to  improve  the  volume  calculation  is  to  remove  a  small 
portion  of  the  charge  near  the  origin.  Th\is,  the  charge  will  be  composed  of  a 
single  8-node  element.  If  the  inside  0.8  in.  of  the  radius  is  removed,  the 
volume  computed  by  DYNA30  is  14  percent  less  than  the  desired  /olume  of 
explosives.  This  is  because  there  is  very  little  volume  inside  of  the  portion 
removed,  and  the  voliuie  is  computed  nuch  more  accurately  for  the  6-node 
element.  Analyses  were  performed  usxng  this  method.  The  deformations  were 
spherically  symmetric,  emd  the  stresses  and  velocities  agreed  reasonably  well 
with  test  data. 
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12.  One  problem  with  using  this  method  is  that  it  requires  fixing  all 
of  the  nodes  on  the  inside  radius  of  the  charge.  This  causes  no  problems  in 
these  analyses.  In  the  SST  study,  there  will  be  a  structure  opposite  the 
charge  on  one  side,  but  there  will  be  no  structure  on  the  other.  This  lack  of 
symmetry  could  mean  that  the  center  of  the  charge  could  move.  If  all  of  the 
inside  nodes  are  fixed,  this  is  not  possible. 

13.  Analyses  were  performed  to  look  at  the  effects  of  using  a 
cylindrically  synunetric  grid  as  compared  to  a  spherically  symmetric  grid.  In 
these  calculations  a  6-node  solid  was  used.  The  volume  of  this  element  is 
computed  correctly  by  DYNA3D,  and  cylindrically  synunetric  deformations  occur. 
Thus,  it  appears  that  the  6-node  element  is  performing  satisfactorily. 
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Fiaure  A1.  Solid  elemente  available  In  DYNA3D 


